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ABSTRACT 

We perform a linear stability analysis of a dust layer in a turbulent gas disk. Youdin 
(2011) investigated the secular gravitational instability of a dust layer using hydrody- 
namic equations with a turbulent diffusion term. We obtain essentially the same result 
independently of Youdin (2011). In the present analysis, we restrict the area of interest 
to small dust particles, while investigating the secular gravitational instability in a more 
rigorous manner. We discuss the time evolution of the dust surface density distribution 
using a stochastic model and derive the advection-diffusion equation. The validity of 
the analysis by Youdin (2011) is confirmed in the strong drag limit. We demonstrate 
quantitatively that the finite thickness of a dust layer weakens the secular gravitational 
instability and that the density-dependent diffusion coefficient changes the growth rate. 
We apply the obtained results to the turbulence driven by the shear instability and 
find that the secular gravitational instability is faster than the radial drift when the gas 
density is three times as large as that in the minimum-mass disk model. If the dust 
particles are larger than chondrules, the secular gravitational instability grows within 
the lifetime of a protoplanetary disk. 

Subject headings: methods: n-body simulations, planets and satellites: formation 

1. Introduction 

Planetesimals are solid objects that are assumed to exist in protoplanetary disks. It is thought 
that they had sizes larger than a kilometer in the solar nebular, which is large enough for gas 
drag forces to be insignificant compared to gravitational scattering. The terrestrial planets and the 
cores of gas giants are considered to have been formed by the collisional accretion of planetesimals. 
However, the planetesimal formation process in protoplanetary disks remains controversial. 



^ Center for Computational Astrophysics, National Astronomical Observatory of Japan, Osawa, Mitaka, Tokyo 
181-8588, Japan 

^ Division of Theoretical Astronomy, National Astronomical Observatory of Japan, Osawa, Mitaka, Tokyo 181- 
8588, Japan 

Department of Physics, Graduate School of Science, Nagoya University, Furo-cho, Chikusa-ku, Nagoya, Aichi 
464-8602, Japan 



- 2 - 



Small du st grains ranging in radii from a few microns to a few centimeters, grow by collisional 
sticking (e.g.. iDominik &: Tielenslll997l : IWurm fc Blumlll998l ). However, the sticking force between 
dust aggregates having radii of larger than one centimeter is so weak that the growth of such 
aggregates may stall. Dust collisions lead to fragmentation of dust aggregates because of the large 
collision velocity (e.g.. iBlum &: Wurmll2000l ). Moreover, meter-sized dust aggregates are strongly 
affected by gas drag and drift radially on account of the hea dwind from the slower rotating gas. 



The timescale of the drift is approximately 100 years at 1 AU (lAdachi et al 
1973) • 



1976 



Weidenschilling 



If the gas flow is laminar, the dust aggregates settle toward the mid-plan e, and the dust 
layer becomes very thin. Such a dense layer may be gravi tationally unstable (jSafronovl Il972l : 
Goldreich &: Ward Il973l : ISekiyal 119831 : lYamoto fc Sekiyal 12004 ) . D ust-rich gas c lumps are formed 



by th e gravitational instability (GI) on the dynamical timescale (ISekiya 



1983 



Yamoto Sekiva 



2006ll. In the clum p, dust aggregates sediment to its center and a planetesimal forms (jSekiva 



198 



i 



Cuzzi et al.ll2008l ) . The planetesimal formation timescale in the clump is l onger than th e dynamical 
timescale, but its size is expected to be the same as the classical estimate (ISekiyalll983l ). Therefore, 
the gravitational instability model has been considered as a promising mechanism of planetesimal 
formation. 

However, if the gas is turbulent, the dust aggregates are stirred up and the dust layer is not 
gravitationally unstable. Various mechanisms for driving t urbulence have b een proposed, includ- 



i ng magneto-rotationa l instabilitv (IBalbus fc Hawlev 



(jWeidenschilling 



1980 



1991 



Sano et al.ll2000l ) and shear instability 



Sekiya h Ishitsul l2000l ) . The resulting turbulence prevents the dust layer 



from settling, and the GI may not occur. 



One possible solution to this problem is the streaming instability model (jYoudin Goodman 



20051 ) . When dust particles are moderately coupled to gas, they clump strongly e n ough to undergo 



gravit ational collapse due to the streaming instability (jJohansen et al.ll2007l . l2009l ). iJohansen et al 
(|2007l ) found that large planetesimals, i.e., Ceres-sized planetesimals, can form directly by the 



streaming instability and subsequent GI. Planetesimal formation by the streaming instability re- 
quires dust aggregates in the protoplanetary disk to become decimeter-sized particles. However, 
the mechanism by which dust particles can grow to decimeter-size is unknown. If dust partic les are 
small and the drag force is strong, streaming turbulence develops (jJohansen &: YoudinI 120071). The 
condit ion for the dust clumping was investigate by the numerical simulations in detail (jBai fc Stone 
2010al lbl). The conditi on depends on the dust size distribution, the dust-to-gas mass ratio and the 



gas pressure gradient (jBai fc Stonell2010bl ). The clumping due to the streaming instability is feasible 
for the small radial gas pressure gradient, the large size dust and the large dust surface density. 



Sekival (jl998l ) investigated the dust density distribution in the turbulence driven by shear 



instability and found that the dust density is larger than the Roche density, the critical density 
for GI, when the dust-to-gas mass ratio is much larger than the solar abu ndance. Thus, i f dus t 
density enhancement mechanisms exist, planetesimals can form due to GI. lYoudin &: Shul (|2002l ) 
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and lYoudin fc Chianei (j2004l ) investigated the evolution of the dust distribution considering the 
radial drift of dust particles. Since the inward mass flux decreases with decreasing radial distance 
from the Sun, dust particles tend to pile up, which may increase the dust-to-gas mass ratio. 



The secular GI can also enhance the dust-to-gas inass ra tio (jWardlll976l : ICoradini et al.lll98ll : 
Wardll200Cl : iGoodman &: Pindoiil200Cll : iMichikoshi et al.ll201Cll ). If we consider the drag force due to 
gas, the dust layer is always secularly unstable, even when Q > 1, where Q is Toomre's Q value 
(jToomrdll964l ). The growth timescale of the secular GI is longer than that of the GI without gas 
drag. Previous studies adopted a simple hydrodynamic model. However, this formulation may not 
be able to correctly handle the effects of the turbulence, and whether the secular GI can develop 
in turbulence remains unclear. If strong turbulence exists, dust particles may be stirred, and the 
dust layer may not collapse. 



YoudinI (j201ll ) extended previous research efforts by considering the turbulent diffusion term. 



He performed a linear stability analysis using an isothermal hydrodynamic equation with a turbulent 
diffusion term. This equation can be applied to a wide range of parameter space. He found that 
the secular GI can develop even when Toomre's Q is greater than unity. When the gas drag force is 
strong, the turbulent diffusion term is signific ant. The added diffusion t erm plays a crucial role. He 
used the a model for the turbulent diffusion (IShakura &: Sunyaevlll973l ). Comparing the timescale 
of the secular GI with the disk lifetime and the radial drift time, he derived the upper limit of a 
value, which is approximately 10~^ to 10~^. Using the turbulence model, he obtained the threshold 
metallicity Z. If dust particles cannot grow large, the secular GI may be more important than the 
streaming instability. 

We have investigat ed the particl e density evolution in turbulence and have derived a similar 
result independently of lYoudinl (j201ll ). In the present study, we re strict our inve stigation to small 



dust particles with a very strong drag force. Although the results of lYoudinI ()201ll ) are more general 



than the results of the present study, we deri ve time evolut ion equations for small dust particles 
in the turbulence in a more rigorous manner. lYoudinl (120111 ) used an intuitive model for the basic 
equations of h is study , which is a hydrodynamic model with a turbulent diffusion term. We verify 
the results of lYoudinl (j201ll ) for small dust particles. If the dust particles are small, the basic 
equation is described by a single advection-diffusion equation. In this case, the dispersion relation 
and eigenfunction can be calcula ,ted in a simple analytic form. Thus, the results of the present 
study are easier to understand. lYoudinl (j201ll ) used a simple turb ulence inodel to calculate the 
threshold metallicity. We adopt the shear turbulence model given by ISekival (j 19981 ) and discuss the 
possibility and condition for the secular GI in shear turbulence. 

The remainder of the present paper is organized as follows. In Section [2l we derive the time 
evolution equation for the dust density distribution for small particles from the Langevin equations. 
This formulation can handle the turbulent diffusion correctly. In Section [3l using the derived time 
evolution equation, we perform a linear stability analysis of the dust layer in turbulence. In Section 
m we describe the shear turbulence model adopted herein. In Section [Sj we calculate the growth 
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timescale of the secular GI and discuss the possibihty of the secular GI. In Section [6l we summarize 
the results. 



2. Time Evolution Equation 



In this section, we discuss the density evolution of small dust particles in a turbulent disk. 
Particles diffuse in phase space because of the random for ce from turbulence. The g eneral formula 
of the diffusion coefficient due to turbulence was given by lYoudin k. Lithwickl (|2007l ) . We consider 
an external force field and diffusion in the radial direction as well as the azimuthal direction, which is 



not considered in lYoudin &: Lithwickl (j2007l ). The external force field does not change the diffusion 
coefficient but determines the drift term if the force gradient is small. The strict derivation of the 
time evolution equation of the surface density with the external force field in the general parameter 
will be investigated in a future study. We assume a thin dust layer and neglect vertical motion of 
dust. 



2.1. Equation of Motion 

We investigate the dynamics of dust particles in turbulence using the stochastic differential 
equation. We consider very small dust aggregates, the stopping time due to gas drag of which is 
very short, ^K, where Tk is the orbital period. 

Using the local Cartesian coordinate system we measure the particle velocity {vx,Vy) 

relative to the local Keplerian shear: Vx = dx/dt, Vy = dy/dt + (3/2)Ox, where Q is the Kepler 
frequency. 

The equations of motion for dust particles are 

dvx 



dt 



2nvy - 'j{vx - Vgx) + fx, (1) 



^ = -^nVx-7{Vy-Vi.y)+ fy, (2) 

where {fx,fy) is the external force that depends on the particle position {x,y), and {vgx,Vgy) is 
the turbulence velocity. The drag coefficient 7 = l/tg is constant. The gas velocities are stochastic 
variables, and Equations ([T]) and ^ are stochastic differential equations and are referred to as 
Langevin equation in statistical physics. 

We assume that the stopping time due to gas drag is shorter than the Kepler time and the 
timescale over which the turbulence velocity changes. Then, the dust velocity is approximated by 
the terminal velocity. The terminal velocity is given by the force balance: 

Vgx 
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2.2. Turbulence Model 

We assume that the ensemble mean of turbulent gas velocity is the Keplerian velocity. Thus, 
we obtain (vgx) = {vgy) = 0, where the angled bracket denotes an ensemble average. 

We consider the turbulence velocity v{t) to be a random time series. The auto-correlation is 
the cross-correlation of the time series with itself (j){t) = {v{to)v{to + t)). If v{t) is stationary, (p{t) 
does not depend on the time Iq. In general, (j){t) is a decreasing function of t. The typical decay 
time of (f){t) is the correlation time tc- If we assume the Kolmogorov spectrum , the auto-correlation 



is described by the following exponential function (jYoudin Lithwicld 120071 ) : 



= a2exp(^-^ ) , (5) 

where a is the velocity dispersion of turbulence. 

We assume that the cor relation time tc is equal to the eddy turnover time te- The typical eddy 



turnover time is ~ ^ ^ (ISekiyalll998l ). We define the gas diffusion coefficients as follows (See 
Appendix A): 

^gyy ~ ^yy^e^ 

^gxy — ^xy^ei (^) 

where Uxx and ayy are the turbulence velocity dispersion, and axy is the magnitude of the cross 
correlation. 



2.3. Advection-DifFusion Equation 

The change in the position of particles due to turbulence is assum ed to be small. Then, the 



time evolution equation for the surface density y, t) is given by (e.g.. lBinnev &: Tremainell2008l ) 
QT, 8 d 1 92 I q2 q2 

where Vx and Vy are the advection velocities, and Dxx^ ^yy^ and Dxy are the diffusion coefficients. 

We assume that the gradient of the external force is small. We obtain the following first-order 
coefficients, which are the drift terms (see Appendix 

^ 1,^ (Ax) _ 7/. + 2n/ 

At^o Ai 72 + 172 ' V ; 
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y At^o At 2(72 + 172) 2 ^ ^ 

where Ax = x — xq and Ay = y — yo are the changes of x and y over the time interval At = t — to, 
to is the initial time, and (xo,yo) is the initial position. We assume that At is smaller than the 
typical timescale At < T but is larger than the eddy turnover time At > tg. These velocities are 
equal to the terminal velocity of dust without turbulent stirring. 

The second-order coefficients are as follows (see Appendix [A|) : 

_ 1 {Ax') _ 7^ZJg.. + 4n^'D,xy + ^n^l'D^yy 

^"""a™0 2 At " (72 + ^2)2 ' ^'^) 



_ 1 {Ax Ay) _ -j^nP^,, + 272(72 + n')D^^y + Anj^D^yy 
^"^"a™0 2 At " 2(72 + 172)2 • 

These coefficients are the diffusion coefficients due t o turbulent stirring. The r adial diffusion coef- 
ficient Dxx is the same as the diffusion coefficient in lYoudin &: Lithwickl (j2007l ). 

If gas drag is strong, then dust particles are well coupled with the gas and the diffusion 
coefficient for the dust is the same as that for the gas 7 Dxx — Dgxx, Dyy ~ D^yy^ Dxy — Dgxy 
We can use this approximation for small dust particles that are much smaller than Im. 

When we consider the external force as the self-gravity of the dust layer, the external forces of 
the X and y components are 

- -I' <^^' 

where is the gravitational potential. The gravitational potential is given by the Poisson equation: 

92 52 92 \ 

7^ + TT^ + 7^ = 4vrGp, (17) 



5x2 Qyl Q^2 



where p is the dust density. 



m 



2.4. Comparison vi^ith Youdin (2011) 

The t i me ev olution equation derived here is the time evolution equation for the special case 



YoudinI (|201ll ). in which the hydro dynamic equation was considered and an additional diffusion 



term was added to the equation of continuity. Since the ga s drag i s assum ed to be strong, the 
Coriolis force is negligible. The time evolution equations in lYoudinl (j201lh without the Coriolis 
force are 

as d{T.Vx) 



dx 



D 



9x2 ' 



(18) 



dt ' s 5x ' 



(19) 



where c is the velocity dispersion, and D is the gas diffusion coefficient. 

If the gas drag force is strong, then the mass flux is estimated in terms of the terminal velocity: 



C2 S/; 



+ 



7 9x 7 

Substituting Equation (pOj) into Equation (fT8|) . we obtain 



as ^ 0(£/./7) 
dt dx 



D, 



eflf 



where -Dgfi is the effective diffusion coefficient: 



D + 



7 



(20) 



(21) 



(22) 



The model of the velocity dispersion and the diffusion coefficient in lYoudinl ( 20 111 ) is as follows: 

(23) 



D 



l + ri 

1 + Ts + At^ 
(1 + t2)2 



(24) 



where Tg = fits, and is the gas diffusion coefficient. For the strong gas drag limit, the velocity 
dispersion and the diffusion coefficient are c ~ (L'gO)^/^ and D ~ Dg. Thus, we obtain 



Deff = D„ 



[I + ^s) ~ I^g. 



(25) 



The effective diffusion coefficient is equal to the diffusion coefficient that we derived for the case in 
which tg <C te- 

If we replace the advection teri n by the terminal velocity in the equation of continuity with a 
diffusion term, then the equation of lYoudinl (j201ll ) is identical to that of the present paper. When 
the gas drag is strong, by solving the equation of motion, the velocity is approximated by the 
terminal velocity with the pressure gradient. Since the pressure gradient term is small compared 
to the diffusion term in the strong drag limit, the terminal velocity is equal to that in the present 
paper. The basic equation of the present paper is derived from the strict stochastic model. Thus, 
the formulation in lYoudinl (|201ll ) is valid at least if the gas drag is strong. 



If the gas drag is weak, the formulation of the present paper, i.e., the advection-diffusion 
equation, breaks down. However, the hydrodynamic equation with a diffusion term for the weak 
drag limit has not yet been verified. When the gas drag is weak, the relaxation of the velocity 
distribution is very slow. In this case, the hydrodynamic picture may break down. The time 
evolution equation for the weak drag limit is arguable. Although we may qualitatively understand 
the physical nature of the instability using the hydrodynamic equation with a diffusion term, in 
order to discuss the physical nature of the instability quantitatively, we should investigate the time 
evolution equation carefully. This issue will be discussed in detail in a future paper. 
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3. Linear Stability Analysis 

3.1. Infinitesimally Thin Disk 

We consider only axisymmetric modes in order that a normal mode analysis can be performed. 
We hereinafter omit the subscript x. From Equations ([9]) and (jlOp . the basic equations are 



V 



if 



72 + 02 



/ 
7' 



(26) 



(27) 



where 7' = 7 + i^^/7 is the effective drag coefficient. For the strong gas drag 7 3> il, we obtain 
7' — 7. 

We perform the linear stability analysis on axisymmetric modes. The unperturbed state is 
uniform in space and time. We consider the perturbation of the form of exp{—i{ujt — kx)), where k 
is the wave number and w is the frequency. The unperturbed and perturbed quantities are denoted 
by the subscripts '0' and '1', respectively. 



Using the thin disk approximation, we can solve the Poisson equation (|Toomrelll964l ). 

27rGSi 



The drift velocity of the x component is 



The linear advection-diffusion equation is 



\k\ 
-ik(j)i 



iwSi = —ikViTjQ — k DT^i. 



(28) 



(29) 



(30) 



Introducing the growth rate fj, = u/i, we obtain the dispersion relation: 

27r|A:|SoG 



-Dk^ + 



7' 



(31) 



If the wave number k is less than the critical wave number 

27rSoG 



kr 



(32) 



the dus t layer is unstable. Th i s feature is sim ilar to the secular GI for the hydrodynamic equation 
model (jMichikoshi et al.ll20ld : lYoudinll201ll ). The dispersion relation indicates that the unstable 
mode always exists, regardless of the diffusion term. The turbulent diffusion makes the growth rate 
small, but it cannot stabilize all modes. The diffusion term is effective for the short wavelength 
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mode. However, the long wavelength mode is still unstable, even though the turbulent diffusion is 
included. 



The maximum growth rate is 



vr 



/^max - j^^i2 ) (33) 



and the wave number of the maximum growth rate is 

"^max — 2 — D'y' ' \ I 

The maximum growth rate is a decreasing function of the diffusio n coeffic i ent. This maximum 
growth rate and the wave number are identical to Equation (56) in lYoudinl (|201ll ) with D = Dg. 



This indicates that the analysis of lYoudinl (|201ll ) is valid when <C t, 



The physical meaning of the secular GI can be understood in terms of the response s to the 



density perturbation. This is the diffusive support secular GI explained in Section 2.2 of lYoudin 



(|201ll ). The responses of the potential and the velocity to the density perturbation are shown in 
Figure [H We assume a sinusoidal density fluctuation. The gravitational potential (pi is induced 
by the density fluctuation. Due to the resultant potential gradient, the matter moves to the local 
maximum point of the surface density at terminal velocity. Therefore, the density fluctuation 
increases monotonically. 

The timescale of the secular GI and the stabilization by the diffusion can be understood as 
follows. First, we neglect the effect of the diffusion. We suppose a wave with wavelength L. 
The speciflc force at a point located at distance L from an infinite straight line with line density 
^ is ~ Gfi/L. We approximate the gravitational field from the wave as the gravitational field 
from an infinite straight line with line density fi ~ LSi. From the force balance between the 
gravitational force and the drag force, the terminal velocity is given as Vi ~ GS1/7'. The resulting 
mass flux is T,Vi ~ GS1S0/7'. Therefore, the timescale of the increase in the surface density 
is Si/(GSiSo/(L7')) — L^' /GTiQ. This corresponds to the growth time calculated by Equation 
(j3ip with D = 0. This mode should be stabilized if the diffusion time L?' /D is shorter than the 
timescale of the secular GI. Thus, the critical length scale for the secular GI is Lcr ~ /GTiq. 
This estimation is consistent with the linear stability analysis. The diffusion can smooth out only 
small-scale density fluctuation. The long wavelength mode is not stabilized by the diffusion. 



3.2. Finite Thickness Disk 

For the sake of simplicity, the dust layer is assumed to be infinitesimally thin. However, in 
reality, the dust particles are stirred by the turbulence, and the dust layer has a finite thickness. 
We account for the effect of the finite thickness of the dust layer by introducing a softening term 




Fig. 1. — Eigenfunctions of Si (solid curve), (dashed curve), and Vi (dotted curve). These 
eigenfunctions are normalized by amplitude. 
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m 



( Vandervoort 


1970: 


Shu 


1984: 



where h is the thickness of the dust layer. This expression indicates that the finite thickness reduces 
the effect of the self-gravity by a factor of 1/(1 + \k\h). 

The dispersion relation for the finite thickness dust layer is 

The critical wave number for the finite-thickness disk is smaller than that for the infinitesimally 
thin disk: 

= Tl^f^^^Y^^- < hr (36) 
The wave number for the most unstable mode fc^ax is given by dfi/dk = 0: 

U2 U h. 

7,13 I o ^max I ^max ^cr _ p. /o^n 
fcmax + ^-^ + -^ ^-U. [6^) 

We calculate the leading term of the power series expansion with respect to h and obtain the 
asymptotic solution of Equation (j37|) . When hkcr ^ 1, ^max approximated as kcr/2. This limit 
corresponds to the infinitesimally thin approximation. When hkcr ^ 1, ^max is approximated 
as kcr:/{2^^^{hkcr:)'^^^). The wave number for the most unstable mode is a decreasing function of 
thickness h. As shown in Figure [21 k'^^^ is smaller than the most unstable wave number for the 
infinitesimally thin disk fc^^x — We can prove that the real positive solution A;^ax/^cr must 

be smaller than 1 /2 for any h. The effect of the finite thickness is to elongate the most unstable 
wavelength. 

The bottom panel of Figure [2] shows the maximum growth rate of the secular GI /i[nax ^ ^ 
function of the dust layer thickness h. When the dust layer is thick, the effect of the reduction 
in the self-gravity due to the dust thickness is significant. Thus, the maximum growth rate is a 
decreasing function of the dust layer thickness. If the dust layer is much thinner than the critical 
wavelength ~ 1/kcr, the effect of the thickness is negligible. Therefore, the thickness factor hkcr 
determines the effect of dust layer thickness. 



3.3. Density-dependent Diffusion Coefficient 



We have neglected the effect of back-reaction of the dust on the gas. If we consider the 
back-reaction, the turbu lence velocity n i a.y cha nge with the dust surface density. According to a 
numerical simulation by iJohansen et al.l (|2009l ) , the particle random velocity in the dense region 
is small, and the resulting diffusion coefficient in the dense region is small. Thus, the diffusion 
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Fig. 2. — Wave number {top panel) and growth rate (bottom panel) for the most unstable mode 
as a functfon of hkcr- In the top panel, the dashed and dotted lines represent kcr/i^^^^ihkcr)'^^^) 
and /ccr/2, which are asymptotic solutions of Equation (I37p . 



- 13 - 



coefficient D may be a decreasing function of the surface density S. In general, D \s a function of 
S. Here, we investigate the density-dependent diffusion coefficient D{Tj). 

The diffusion coefficient is assumed to depend on the surface density: 

L> = D(S)~D(So + Si)~ A + (38) 
where Dq = D{T,q). The diffusion term in the advection-diffusion equation is 

^ ^-'^o + :HtSoU-^. (39) 



dx^ \ dT, J dx^ 

In other words, the term obtained from dD/dT, is added to the diffusion coefficient. Thus, when 
we consider the density-dependent diffusion coefficient, we must replace the diffusion coefficient in 
Equation ([5T]) by 

dD ( d\ogD\ , , 



Therefore, we obtain 



where 



, = _Z,„(l_^,,= + 5dAE!!^, (41) 

7 

/3 = -^- (42) 



dlogS 

If we adopt the power-law model, we obtain D cx 

When /? > 1, any mode is unstable, although self-gravity is not considered. Essentially, this is 
the same as the diff usion instability (viscous instability) that was discussed in explaining the ring 



structure of Saturn ([Lukkaril Il98ll : iLin &: Bodenheimeii llDBll ) . If /3 > 1 , SZ) is a local minimum 



value when the density is a local maximum value. Thus, the density at the local maximum point 
of the surface density increases due to the resultant diffusion flux. When /3 = 1, the diffusion term 
vanishes. In this case, due to self-gravity, any mode is unstable. When /? < 1, the diffusion term 
stabilizes the short wavelength modes. Thus, the maximum growth rate exists. The maximum 
growth rate and the wave number are 

_ 7r2£2G2 

(l_;3)Z)oy2' ^ 

Do(l-/3)y ^ ^ 

The growth rate increases with /3 because the effective diffusion coefficient decreases with /3. 

In the present paper, we focus on the shear turbulence in the application. The diffusion 
coefficient due to shear turbulence does not depend on S. Thus, in the following, we will not 
discuss the secular GI for the density-dependent diffusion coefficient. Note that this effect can be 
important in some types of turbulence. In particular, when /? ~ 1, the maximum growth rate 
depends sensitively on /?. 
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4. Shear Turbulence Model 

The hnear stabihty analysis reveals that a small density fluctuation can grow in a turbulent 
disk regardless of the turbulence strength. We apply the secular Gl to the turbulence due to the 
shear instability in a protoplanetary disk. 



4.1. Minimum-mass Disk Model 



The surface d ensity is estimated based on the minimum-mass disk model (jHavashi 
Hayashi et al.lll985l ). The surface densities of gas and dust are 



1981 



a \-3/2 



g/cm^, 



Sd = 7.1/d(^) ^"g/cm2, 



(45) 



(46) 



where a is the distance from the Sun, and /d and /g are enhancement factors of the surface densities 
of gas and dust. In the minimum-mass disk model, we assume /d//g = 4.2 outside 2.7 AU because 
H2O condenses. The gas temperature T, the gas density pg on the mid-plane, the sound velocity 
Cg, and the pressure gradient r) are obtained as follows: 



Mm) 



-1/2 _ 



I.4XI0-V.(^)-"./ 



1.1 X 10^ 



lAU 



K 

11/4 
-1/4 



cm 



cm/s. 



a \i/2 



1.81 X 10-3 (^) 



(47) 

(48) 
(49) 
(50) 



4.2. Structure of the Dust Layer 



The dust particles are assumed to be very small. In this case, the gas drag force is strong 



As dust settles toward the niid-plane, shear instability is inevitable (ISekiya &: Ishitsul I2OO01 . 12001 



Ishitsu Sekivall2002l . l2003l : lMichikoshi &: Inutsukall2006l ). The resultant turbulence prevents dus t 
from settling further. Consequently the dust layer attains a quasi-equilibrium state (|Sekivalll998l ). 



The stability of the shear flow is determined by the Richardson number J: 



J 



g dp/dz 
p {dv/dz)' 



(51) 
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where g is the vertical gravitational acceleration. If J is less than the critical value Jcr, the shear 
flow becomes turbulen t. The linear stability analysis reveals that the critical value is Jcr = 0.25 
( Chandrasekhaii Il96lh . Numerical simulation without the Coriolis force confirmed the critical 



Richardson number to be approximately 0.25 (lBarraiicdl2009l ) . Numerical simulations that include 
the effect of rotation indicate that Jcr — 0.1 (jchiang 2008 ). In the case of small dust, the critical 
Richardson numb er is Jcr — 0.2 for the solar metallicity disk, which increases with metallicity 
(jLee et al.l l2010bl ). Therefore, the critical Richardson number is Jcr = 0.1 — 0.25. When we 
estimate various quantities, we use Jcr = 0.1, which is an optimistic assumption. 



Sekival (|l998l ) argued that the dust density distribution is characterized by assuming that 



the Richardson number J{z) is equal to the critical R ichardson number Jcr everywher e. This 



assumption was explored by the numerical simulations (jBarranco 



2009; Lee et al. 



2010 



30). The 



Richardson number of the dust layer in a protoplanetary disk is (jSekiyalll998l ). 



J{z) 



z + 



JO 

where pd{z) is the dust density. Solving J{z) 



{pg + Pd{z))dz 



{mPgfdPd/dz' 



(52) 



Jcr, we can calculate the dust density distribution. 



The thickness of the dust layer is (jSekivalll998l ) 

/id = 



Jcr^/a\/l - E?, 

where R is the ratio of the gas density to the total density at the mid-plane: 



R 



Pg + Pd{^) 



The ratio R is obtained as the solution of the following equation (jSekivalll998l ): 

S 



l + q 



log 



l + v/l-((i? + g)/(l+g))^ 
{R + q)/{l + q) 



^l-{{R + q)/{l+q)Y, 



(53) 



(54) 



(55) 



where S = Sd/(2-v/ Jcr??«Pg), and q = 47rGpg/r2^ . The parameters S and q in the minimum-mass 
disk model are 

5 = 0.030/.(-) 



s 



(56) 
(57) 



When S = 0.3 and q = 0.03, we obtain R = 0.60 , and thus Pd/Pg — 0.68. Since the dust 
density is mu ch smaller than the Roche density (|Sekiyalll983l ). the dynamical GI does not occur. 
Sekival (|l998l ) reported that the dust density can be larger than the Roche density if the dust-to- 
gas mass ratio is much larger than the solar abundance, fd/fg ^ 20, for T = 280 K at lAU. The 



critical dust-to-gas mass ratio depends on the disk parameters (jYoudin &: Shull2002l ). If the disk 
temperature is 100 — 170 K at lAU, the critical dust-to-gas mass ratio is about 2 — 10. In the 
present paper, we assume that the dust-to-gas mass ratio is the solar abundance fd/fg ~ 1 and 
investigate the possibility of the secular GI. 
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4.3. Diffusion Coefficient 



Setting the mass flux due to the sedimentation velocity equal to the diffusion flux, we have 



( Cuzzi et al 



1993 



Sekivalll998l ). 



Ddpd 



Pd dz 

From Equations (j52p and (jSSp . we obtain the diffusion coefficient D, which depends on z, 

plPdiz) 



D{z) = Jcrts-rfvl^ 



where F„ is the effect of the self gravity and has the following upper limit: 



1 + ^ 
z .10 



l + ^]dz<l + . 

pg. 



1 + 



MO) 

Pg 



1 + -. 
R 



(58) 



(59) 



(60) 



Since q/R <^ 1 is satisfied in the standard model, we assume that Fsg — 1- Accordingly, the 
diffusion coefficient is 

^ 7.2 2 plPdjz) 
D{z) ~ JcAf? ^ 61 

If Pd(0)/pg > 1/2, D{z) has a maximum value when p^{z)/ pg = 1/2, then: 



We assume the typical diffusion coefficient in the dust layer to be equal to half of D„ 

2 



(62) 



(63) 



This is the vertical diffusion coefficient of the turbulence. The vertical diffusion coefficient is 
assumed to be equal to the radial diffusion coefficient. 



YoudinI (1201 ih used D ~ tsV^V^ (Equation (68) inlYoudinI (120111 1). The diffusion model used 



herein has the same dependence as the model of lYoudinl (|201ll ). except for Jcr, and the coefficient 
of the diffusion model used herein is much smaller. 



5. Secular Gravitational Instability of a Dust Layer 

5.1. Timescale and Wavelength 

In this section, the growth time of the secular GI in a protoplanetary disk with shear turbulence 
is estimated. First, we examine the dust layer thickness. We can calculate the thickness factor from 
Equation ([53]) as follows: 

hkcr = v^W^cr Vl-R^ < y/j^Tjakcr = 0.12/d ( ^ j ■ (64) 
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The thickness factor does not depend on the dust particle radius and the distance from the sun 
because they are canceled out. Since the thickness factor is less than unity in the standard model, 
we adopt the infinitesimally thin disk approximation in calculating the growth time of the secular 
GI. 



The growth time of the secular GI in the thin disk approximation is 



IT- 



(65) 



There are two models of gas drag: the Stokes drag and the Epstein drag. If th e dust particle 



radius is less than the mean free path of gas, the Epstein drag is appropriate (e.g., lEpsteinlll924l : 



path, t he Stokes drag is a' 



IS e.E 



Adachi et al. 



197' 



propriate (e.g.. 


Stokes 


1851 


Nakaeawa et al. 


1986|) 



11/4 



cm. 



Thus, when a < Ocr , the Stokes drag should be used, where 

acr = ft''' 



Adachi et al.lll976l ). The mean free path I 

(66) 



1cm 



and Ocr is determined by Z = rp. The stopping time is 



(67) 



a > Or 



where pp is the bulk density of the dust particles, and is the radius of the dust. The growth 
time of the secular GI is 



tGl 



T ( :h 

/g ( Ja 
'd 



2.4 X 10 



1.6 X 10^41 I — 



n VO.l 



Pp ^ 




3g/cm'^y 


V 1mm 




Pp ^ 




3g/cm'^y 


\ 1mm 



a \ii/4 



lAU 



-1 



years 



years {a < 

(a > Oc 



(69) 



The growth time for a > Ocr does not depend on the distance from the Sun a. 
The critical wavelength of the secular GI is 



An 



27r 



Di _ 2Jcrr?^4 _ q n V in-2 f-i ^ ^" 



/ a \3/2 



0.1/ VAU/ 



AU. 



(70) 



This does not depend on the radius of the dust particle. 
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5.2. Instability Conditions 



5.2.1. Comparison with the Radial Drift 

If the radial drift of dust particles due to gas drag is faster than the secular GI, the secular GI 
is inefficient. In stead, the particle p i leup mechanism due to the radial drift enhances the surface 
density of dust (jYoudin k. Shull2002l : lYoudin k. Chianel 12004 ). Here, we compare the timescales of 
the secular GI and the radial drift. 



Youdinl (|201ll ) examined the same condition. Comparing the timescale of the secular GI with 
the ra dial drift, he derived the critical turbulence alpha value Omax (Equation (53) in lYoudin 
(|201ll )). He applied amax to the turbulence model D = tsrfv\ and calculated the metallicity 
thresholds. In this section, we perform es sentially the same analysis. However, we assume the 
turbulence model reported bv lSekival (|l998l ). 



According to the linea r stability analysis, the streaming instability is faster than the radial drift 
(jYoudin Goodmaiill2005l ) . In the present paper, we consider small dust particles. The wavelength 
of the streaming instability is very short for small dust particles. The numerical simulation of tightly 
coupled dust aggregates suggest tha t the concentration caused by the streaming instability is not 
very effective (jJohansen et al.l 120071 ). The clumps are small and short-lived. Thus, in the case of 
small dust particles, the streaming instability contributes to the turbulence, but may not promote 
planetesimal formation. In the present paper, we ignore the streaming instability and focus on the 
radial drift. 



The radial drift velocity of dust particles is (jNakagawa et al.lll986l ) 

If the gas drag is strong Rt^^ ^ 1, the radial velocity is approximated as 

Vr ^ 2R^ts^r]VK- 



The radial drift time of dust particles is 



1 



Vr 2l]R^nHs' 

Thus, the ratio of the timescale of the secular GI to the radial drift is 



277r2G2S 



(71) 



(72) 



(73) 



(74) 



This ratio does not depend on the dust particle radius. As the dust surface density increases or the 
pressure gradient rj decreases, the secular GI becomes faster. The ratio R is positive and smaller 
than unity. 
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We examine the temperature dependence. Since ry and R depend on T, the ratio tci/tr depends 
on T and S^. Figure [3] shows the condition for the secular GI in the T-Sj plane at 1 AU. We 
numerically calculate R from Eq. ()55p . In the case of the standard model, the radial drift is 
faster than the secular GI. If the surface density is large or the temperature is low, the secular 
GI can occur. In the case of /d = 1, for the secular GI, the temperature must be lower than 220 
K. As the temperature decreases, the pressure gradient r] decreases. Thus, the vertical shear and 
the resultant shear instability weaken. In the case of T = 280K, if /d > 1.7, the secular GI is 
faster than the radial drift. In general, the realistic surface density can be larger than that of the 
minimum-mass disk model. We investigate the case in which the surface density is larger than that 
of the minimum-mass disk model, i.e., /g, /d > 1- 

We calculate the condition in which the secular GI is faster than the radial drift in the /d~/g 
plane. We define the critical value /d,cr = fd at which tci/tr = 1. When /d > fd,cr{fg, Jcr,a), the 
secular GI is faster than the radial drift. Figure H] shows that /d,cr is an increasing function of /g 
at a = lAU, where we adopt Jcr = 0.1. If /g < 2.2, /d,cr for 1 AU is larger than /g. Thus, the 
enhancement of the dust-to-gas mass ratio is necessary, so that the secular GI is faster than the 
radial drift for /g < 2.2 at 1 AU. On the other hand, for /g > 2.2, although the dust-to-gas mass 
ratio is the standard value fd/ fg = 1) the secular GI can be faster than the radial drift. Therefore, 
if the protoplanetary disk is more massive than the minimum-mass disk, the secular GI can be 
effective at 1 AU, although the dust-to-gas mass ratio is the standard value. Since H2O condenses 
at 5 AU, the standard dust-to-gas mass ratio is 4.2. Thus, /d is larger than /d,cr) although /g = 1.0. 

Next, assuming the standard dust-to-gas mass ratio, we examine the condition in which the 
secular GI is faster than the radial drift in the a-/g plane. We define the critical value fg^criJcr, cl) 
as /d,cr(/g,cr,'/cr,a)//g,cr = 1 for o < 2.7AU and /d,cr(/g,cr, -/cr, a)//g,cr = 4.2 for a > 2.7AU. If 
fg > /g,cr) the secular GI is faster than the radial drift with the standard dust-to-gas mass ratio. 
Figure [5] shows the critical /g^cr as a function of the distance from the Sun a for Jcr = 0.1 and 
0.25. The critical value /g^cr has a maximum value at a = 2.7AU. Therefore, for the case in which 
Jcr = 0.1, if fg > 3, the secular GI is faster than the radial drift in the entire disk. When = 0.25, 
the turbulent diffusion is stronger than that for Jcr = 0.1. Then, a more massive disk, such that 
fg > 5.5, is necessary for the secular GI. 



5.2.2. Comparison with the Disk Lifetime 

The secular GI should occur within the disk lifetime. The disk lifetime is typically 10^ — 10^ 
years. Figure [6] shows the timescale of the secular GI. We adopt /g = 3 and Jcr = 0.1. In the 
inner disk, the Stokes drag model is appropriate because the mean free path of the gas is less than 
the dust particle radius. The timescale of the secular GI is proportional to the inner disk. 

On the other hand, in the outer disk, the Epstein drag model is appropriate because the mean 
free path of the gas is longer than the dust particle radius. The timescale of the secular GI does 
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100 150 200 250 300 

T(K) 

Fig. 3. — Critical surface density for the secular GI at 1 AU as a function of the temperature T 
for /g = 1. If the surface density is larger than the critical value, the secular GI is faster than the 
radial drift. The plus symbol denotes the standard minimum-mass disk model. 
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Fig. 4. — Condition for the secular GI in the /g-/d plane. The solid curve denotes the critical /d,cr 
value as a function of /g for 1 AU. When /d > /d,cr, the secular GI is faster than the radial drift. 

The dashed curve denotes the critical /d,cr value for 5 AU. The dotted line denotes the case in 
which /d = /g, which is the standard dust-to-gas mass ratio inside the snow line. The dot-dashed 
line denotes the case in which /d = 4.2/g, which is the standard dust-to-gas mass ratio outside the 
snow line. 
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Fig. 5. — Critical gas enhancement factor /g^cr as a function of the distance from the Sun a for 
Jcr = 0.1 (soUd curve) and that for Jcr = 0.25 (dashed curve). If /g is larger than the critical 
/g,cr, the secular GI is faster for the standard dust-to-gas mass ratio. For a < 2.7AU, we calculate 
/g,cr from /d,cr//g,cr = 1- For a > 2.7AU, we calculate /g,cr from /d,cr//g,cr = 4.2 because H2O 
condenses. 
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not depend on the distance from the Sun a. The discontinuity at a = 2.7AU corresponds to the 
snow hne. Since the surface density of dust is large outside the snow hne, the secular GI is rapid. 
Since the drag force from gas is weaker for larger dust particles, the timescale of the secular GI is 
faster for larger dust particles. The typical chondrule radius is 1mm. Thus, the growth time of the 
secular GI is 5.3 x 10^ years. The secular GI can grow within the disk lifetime. If dust aggregates 
grow to 1 cm, the growth time is 5.3 x 10^ years, which is very rapid compared to the disk lifetime. 

In Figure[6l the growth time has a maximum value at the snow line a = 2.7 AU. From Equation 
(j671) . the Epstein drag law is appropriate at a = 2.7 AU for 

rp < 15.4/g-^cm. (75) 

Since we assume small dust particles, such as 1 mm, this condition is satisfied. Therefore, the 
maximum growth time of the secular GI in the protoplanetary disk is evaluated in terms of Epstein 
drag. 

The condition whereby the secular GI grows within the disk life time tjisk is given as follows: 

r„ > r„„ = 0..e| [f^) (^) " (^) " mm. (76) 

When /g = /d = 3, the critical radius in order for the secular GI to grow within lO'^ years is 
Tp ~ 0.053 mm. 



6. Summary and Discussion 



We investigated the secular GI for the case in which du s t part icles are very small. A general 
analysis of the secular GI has been performed by lYoudinl (|201ll ). He used the hydrodynamic 
equation with the additional diffusion term. His result may be applicable for even large dust 
particles. However, his formulation is somewhat intuitive. We herein restrict ourselves to small 
dust particles and investigate the same instability rigorously. We obtained esse n tially th e same 



result, although the present study was conducted independently of lYoudinl (j201ll ). lYoudinl (120111 ) 



considered the gerieral t urbulence, whereas we focused on the specific shear turbulence model 
reported by ISekiyal (|l998l ) and discuss the condition for the secular GI in detail. 



In the case of small particles, the stopping time is shorter than the turbulence auto-correlation 
time, ts <^ tc- We can neglect the particle inertia because the gas drag relaxation is very rapid. 
We derived an advection-diffusion equation for small particles in turbulence from the Langevin 
equation considering the external force field. We confirmed that the radial diffusion coefficient 



derived in the present paper is the same as that derived in lYoudin &: Lithwickl (|2007l ). 



Considering the self-gravity of the particles with a thin disk approximation, we performed 
linear stability analysis of the advection-diffusion equation. We calculated the growth rate of the 
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Fig. 6. — Timescales of the secular GI for rp = 1cm (solid curve), rp = 1mm (dashed curve) and 
Tp = 0.1mm (short dashed curve). We assume that /g = 3 and Jcr = 0.1. 
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secular GI. Although the diffusion can stabilize the short wavelength mode, the long wavelength 
mode is always unstable, even though the turbulence is strong. Next, we considered the effect of the 
finite thickness of the dust layer. The thickness of the dust layer weakens the secular GI because 
of the lower density of a dust layer of finite thickness. If the diffusion coefficient depends on the 
surface density, the growth rate changes. In particular, if the diffusion coefficient decreases with 
the surface density, the growth rate of the secular GI can become very large. 

Assuming the minimum-mass disk model, we calculated the growth rate of the secular GI. We 
adopted the constant Richardson number density distribution model as the dust layer model (jSekiva 
19981 ). We calculated the turbulent diffusion coefficient from the density distribution assuming 
the mass flux balance and confirmed that the thin disk approximation is valid when the critical 
wavelength is smaller than the thickness of the dust layer. If the dust-to-gas mass ratio is large or 
the temperature at the mid-plane is lower, the secular GI is faster than the radial drift due to the 
gas drag. In the case of the minimum-mass disk model, the radial drift is faster than the secular 
GI. The realistic surface density can be larger than that of the minimum-mass disk model. If the 
surface density is more than three times as large as that in the minimum-mass disk model, then 
the secular GI is faster than the radial drift. In this case, the secular GI with millimeter-sized dust 
can develop within 5.3 x 10^ years. 

Note that the growth time is of the e-folding time scale, which is the time in which the 
perturbation quantities increase by a factor of e. Generally speaking, the planetesimal formation 
time should be larger than the growth time of the secular instability. The realistic planetesimal 
formation time depends on the amplitude of the initial perturbation and the nonlinear process that 
occurs after the linear growth. Further investigation is necessary in order to clarify planetesimal 
formation by secular GI. 

We restricte d the present study to the strong drag case tg <C te and confirmed that the analysis 
of lYoudinI (120111 ) is valid when <^ te- His analysis may be applicable to more general parameters, 
such as ts ^» te- Strictly speaking, we should solve the time evolution equation of the velocity 
distribution in this parameter regime, because, in general, the relaxation process of the velocity 
distribution is not negligible. In addition, visc osity nia.y not be negligible. The hydrodynamic 
equation with a turbulent diffusion term used bv lYoudinI (|201ll ) should be more rigorously justified. 
We will investigate the time evolution equation and the secular GI for more general cases in a future 
study. 



In the linear stability analysis, we assumed the axisymmetric modes. In general, non-axisymmetric 
mo des can be excited . In fact, non-axisymrnetric density patterns were observed in previous stud - 



ies (jTanga et al 



2004 



Michikoshi et al. 



2007 



Wakita Sekivall2008l : Michikoshi et al 



200S 



2010|). 



These structu res were i nterpreted as being similar to self-gravit ating wakes induced by GI and 



Kepler shear (|Saldll995l ). When we consider the gas-fre e model (jMichikoshi et al. 



2007 



20091 ) or 



when the gas drag force is weak ([Michikoshi et al.ll2010l ). the timescale of the GI is on the order 
of Kepler time. The timescale of the secular GI under strong gas drag and turbulent diffusion is 



- 26 - 



much longer than Kepler time. Thus, the Kepler shear becomes significant relative to the GI. The 
shear makes a non-axisymmetric mode a tight trailing mode, which may eventually be smoothed 
by the turbulent diffusion. Thus, the non-axisymmetric mode may not be able to grow sufficiently. 
The axisymmetric mode is expected to grow in the secular timescale. 

Axisymmetric high-density rings form due to the secular GI. As the secular GI develops, the 
dust-to-gas mass ratio increases, and the dust density at the mid-plane increases. Finally, the dust 
density at t he mid-plane exceeds the Roche de nsity when the dust-to-gas mass ratio becomes larger 



than 2 — 20 (jSekivalll998l : lYoudin &: Shull2002l ). When the dust surface density in the ring increases 



by a factor of F, the ring width decreases, as follows: 



Aa ^ %i = 6.0 X 10^3 (^)(^) AU, (77) 



F -"^ \10 J VO-1 / ^AU/ 

where Amax = 27r/A;max- On the other hand, the critical wavelength of the dynamical GI of the thin 
rotating disk is 



^O. = ^ = 3.2xl0-V.(^5)(^) AU. 



(78) 



The condition for Aa ^ Aqi is written as 



'J a 



1/2 



UF « 43.5 (^^J . (79) 

As shown in Figure El the critical /g value for Jq- = 0.1 at 1 AU is approximately 2. Then, 
the critical wavelength of the dynamical GI is slightly smaller than the shrunken width of the 
axisymmetric dense ring. We may be able to apply the dynamical GI theory of the rotating disk. 
In this case, the planetesimal mass and radius are estimated by the classical GI theory: 

.Hp, ~ EjA^, = 1.63 X lO^VI f^) (^)'"8. 

The planetesimal size is larger than the classical planetesimal because the dust surface density is 
large. 

In general, the factor F may be larger than approximately 10 and the initial value of /d may 
be larger than 2. Then, the classical GI theory of the thin rotating disk breaks down because of 
Agi > Aa. In this case, stability analysis of a thin ring or filament of dust with turbulent stirring 
is necessary in order to estimate the planetesimal size. We intend to examine the stability of a thin 
filament of dust with turbulent stirring in a future study. 



We thank Minora Sekiya for helpful comments on this research. 
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A. Derivation of Coefficients 

We assume the timescale T of the phenomenon consider herein to be much longer than the 
eddy turnover time te, i.e., T ^ tc- In the following discussion, the auto-correlation function is used 
in the integrals. Since the auto-correlation is negligible for t > te, the auto-correlation function 
described by Equation ^ in the integrals is replaced by the delta function: 

m ^ cm, (Ai) 

where 6{t) is the Dirac delta function. Assuming that the integral over the time t is invariant, 

C5{t)dt = j vleicp{-j-^dt, (A2) 

we determine the coefficient of the delta function C as follows: 

C = 2tcvl (A3) 

The auto-correlation functions of the x and y components of the turbulence velocity are 

^^^{t) ~ 2tealj{t) = 2Dg^J{t), (A4) 

(l)yy{t) ~ 2te(jly5{t) = 2Dgyy6{t), (Ab) 

where itx and tty are the turbulence velocity dispersion, Dgxx and Dgyy are the gas diffusion 
coefficients, where Dg^x = teC^xx ^gyy ~ ^e^^yy "^^^ cross-correlation function of the x and y 
components of the turbulent velocity is 

(l)xy{t) ~ 2te(jl 6{t) = 2DgxyS{t), (A6) 



This preprint was prepared with the AAS IATJhjX macros v5.2. 
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where axy is the magnitude of the cross correlation, and Dg^y = ieC^xy 
The formal solution of Equation ([3]) is 

(A7) 

The external forces fx and fy are assumed to be smooth functions with respect to x and y, respec- 
tively. The integrals of the external forces are approximated as follows: 



/ fxdh ^ UAt + ^ / - x{to))dh, (A8) 



where we assume the change in position to be small. If the gradient of the external force is small, 
the second term is negligible. The integral of the external force is assumed to be as follows: 

rto+At 

/ fxdh ^ fx^t. (A9) 

•Jto 



Therefore, the first-order solution is written as 

„2 rto+At 27^1 rto+At 

In the same way, the formal solution of Equation @ is 



2 rto+At f^o+At J , , 

72 02 y, ^ 72 02 72 + J^2 ^2 _^ f^2 



Aw ~ / -^dt = [vy Ox dt 

to dt 4 V' 2 ' 

'"s^dh + -o—^ / ^§2/^*1 + .,2 , 02 ^^ - ^T7:::2— Ff2T^* 



2(72 + ^^2);^^ -8-— ■ ^2 + f^2y^^ -8.-- ■ ^2 + f^2 ^ 2(72+02) 

3 3 720 37O2 f^O+At rt 



2^"°^*-2 7To^i ^'^'^^ " TTo^ lo ^'^l/'^'^^y ^^''^ 

Since the ensemble means of Vgx and Ug^^ are zero, the integrals are also zero: 

rto+At rto+At 
{ / Vgxdti) = / {vgx)dti = 0, (A12) 
Jto J to 

rto+At rto+At 

( / Vgydtl) = / {Vgy)dh = 0. (A13) 

Jto Jto 

Thus, the ensemble means of Equations (lAT(i|) and (fATT]l are 
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We readily obtain the first-order coefficients as shown by Equations (jlOp and ([TT 
The ensemble mean of Ax^ is 

4 i-to+At i-to+At 

(^^') - (^2 ^ ^2)2 / / '^t2(^gx(tl)^^gx.(t2)) 

2^3j^ rto+At i-tQ+At 



2o2 I-to+At I-to+At 

/ (iti / (it2(?^gy(tl)^^gy(t2)>. (A16) 

J to J to 



(72 + 1^2)2 



Since we adopt the white noise approximation, the auto-correlation functions are given by the Dirac 
delta functions as Equations (jA4j) , (|A5P , and ()A6P . The integral in the first term of Equation ()A16P 
is evaluated as 



I-to+At I-to+At I-to+At I-to+At 

/ dh dt2{v^x{ti)v^x{t2)) = dti dt22DgJ{ti-t2) = 2Dg^At. (A17) 

J to J to ^to 'J to 

Thus, we obtain the diffusion coefficients D^x as Equation (I12p . 
The ensemble mean of Ay2 ig 



2^2 I-to+At I-to+At 

i^y') - 4(^2 ^ ^2)2 dt2{v^x{tl)Vgx{t2)) 



j-to+At I-to+At 

dti / dt2{Vgx{tl)Vgyit2)) 



(72 + 02)2 



to ■J to 



2^2 I-to+At I-to+At 
+ (^2^^2)2 / / '^i2(^gy(il)t'gs,(i2)) 



to ■J to 



where 



+A. (A18) 



3 72ri I-to+At i-ti ^^q2 I-to+At rti 

A = {2[ ---^ J^^^ dt, j^^^ dt2V,x - dt, j^^ dt2V,y 

I-to+At ^2 pta+At 



Vgxdti + I Vgydti 



^ 2 



/ 3 Y^ I-to+At i-ti 3^j^2 I-to+At i-ti 

The expansion of A contains multiple integrals, as follows: 

I-to+At I-to+At i-ti 2 

/ dh / dh / dt26{t2 - ts) = -Af, (A20) 

Jto Jto Jto ^ 
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/•to+At i-ti i-to+At i-ts 2 

/ dti / dt2 / dts / dU6{t2 - U) = -At^. (A21) 

Jto Jto Jto Jto 

These integrals are higher-order terms. Accordingly, the term A is ~ O(At^), which is negligible. 
The diffusion coefficients Dyy are then obtained as shown in Equation (|13p . 



